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Abstract 

Climate exhibits a vast range of dissipative structures. Some have characteristic times of a few 
days; others evolve on thousands of years. All these structures are interdependent; in other words, 
they communicate. It is often considered that the only way to cope with climate complexity is to 
integrate the equations of atmospheric and oceanic motion with the finer possible mesh. Is this the 
sole strategy? Aren't we missing another characteristic of the climate system; its ability to destroy 
and generate information at the macroscopic scale? Paleoclimatologists consider that much of this 
information is present in palaeoclimate archives. It is therefore natural to build climate models such 
as to get the most of these archives. The strategy proposed here is based on Bayesian statistics and 
low-order non-linear dynamical systems, in a modelling approach that explicitly includes the effects of 
uncertainties. Its practical interest is illustrated through the problem of the timing of the next great 
glaciation. Is glacial inception overdue, or do we need to wait for another 50,000 years before ice caps 
grow again? Our results indicate a glaciation inception in 50,000 years. 
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L'analyse mathematique pent deduire des 
phenomenes generaux et simples V expression 
des lois de la nature; mais I'application speciale 
de ces lois a des effets tres-composes exige une 
longue suite d' observations exactes. 

Joseph Fourier (1768 - 1830) 

1 In troduction 

This quote by Joseph Fourier appeared first 
in the "discours prehminaire" of the analyti- 
cal theory of heat [1 . At a time when the 
reversible Newtonian equations were champi- 
oned by Pierre-Simon Laplace (1749 - 1827) and 
Joseph Louis Lagrange (1736 - 1813), the ir- 
reversible equations governing heat propagation 
constituted a genuine mental revolution. With 
this sentence, Fourier arguably sets the founda- 
tions of complex system theory. He repeated 
it at least once, to conclude his memoire sur 
les temperatures du globe terrestre et des es- 
paces planetaires [2^ in which Fourier formulates 
what is known today as the "greenhouse effect" . 
Fourier confesses that "the question of Earth's 
temperature is one of the most important and 
difficult of all the Natural Philosophy" [2j and 
solving it was one central motivation for the 
theory of heat. Clearly, Fourier had fully per- 
ceived the complex character of the climate sys- 
tem. How, two centuries later, do we cope with 
climate's complexity? Which mathematical anal- 
ysis is the most appropriate to get the best out of 
observations? With this paper we would like to 
convince the reader that the most complex model 
is not necessarily the most useful. Predicting and 
understanding the climate system requires a con- 
sistency between the level of complexity of of ob- 
servations, model prediction and what one wants 
to predict. Choosing the right model is thus also 

^ Mathematical analysis allows you to deduce nature's 
laws from general and simple phenomena; but applying 
these laws to highly composite effects requires a long series 
of exact observations 



a question of information theory. 

The case will be illustrated through a polemic 
currently taking place in the circle of Quaternary 
climate scientists. Here is it. As we shall see in 
more detail, the climate history of the past few 
million years is characterised by repeated tran- 
sitions between "cold" (glacial) and warm (in- 
terglacial) climates. The first modern men were 
hunting mammoth during the last glacial era. 
This era culminated around 19,000 years ago [3] 
and then declined rapildy. By 9,000 years ago cli- 
mate was close to the modern one. The current 
interglacial, called the Holocene, has lasted long 
enough compared to previous interglacials. The 
polemic is about when it is supposed give way to 
a new glacial inception, keeping aside human ac- 
tivities that have most probably perturbed nat- 
ural cycles. 

On the one side. Professor of Environmen- 
tal Sciences Bill Ruddiman carefully inspected 
and compared palaeo-environmental information 
about the different interglacial periods. This 
comparison exercise let him to conclude that 
glacial inception is largely overdue [U |5]. Ac- 
cording to him, the Holocene was not supposed 
to be that long, but the natural glacial inception 
process was stopped by an anthropogenic per- 
turbation that began as early as 6,000 years ago 
(rice plantations and land management by an- 
tique civilisations). On the other side. Professor 
Andre Berger and colleagues developed a mathe- 
matical model of the climate system, rated today 
as a "model of intermediate complexity" [6l [7] in- 
cluding 15,000 lines of FORTRAN code to solve 
the dynamics of the atmosphere and ice sheets 
on a spatial grid of 19 x 5 elements, with a 
reasonably extensive treatment of the shortwave 
and longwave radiative transfers in the atmo- 
sphere. Simulations with this model led Berger 
and Loutre conclude that glacial inception is not 
due before 50,000 years as long as the CO2 atmo- 
spheric concentration stays above 220 ppmv. [8] 
Who is right? Both (Crucifix and Berger argued 
that the two statements are not strictly incom- 
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patible [S])? None? Both Ruddiman and Berger 
judge that it is possible to predict dimate thou- 
sands of years ahead but is it a reahstic expecta- 
tion after ah? Michael Ghil wondered "what can 
we predict beyond one week, for how long and 
by what methods?" in a paper entitled "Hilbert's 
problem of the geosciences in the XXIst century" 
|10j . This is the fundamental motiviation behing 
the present article. 

2 Steps towards a dynamical 
model of palaeoclimates 

2.1 Some general remarks about com- 
plex system modelling 

A system as complex as climate is organised at 
different levels : clouds, cloud systems, synop- 
tic waves, planetary waves, pluri-annual oscil- 
lations such as El-Niho, glacial-interglacial cy- 
cles. ... It is not our purpose to explain here 
how, in general, patterns emerge in complex sys- 
tems ([m [12] are up-to-date references on the 
subject) but it is useful to have a few notions in 
mind. Complex systems and their components 
act as information processors. This means that 
their dynamics is such that they can destroy, am- 
plify and even create information. The difficult 
mental barrier to overcome for physicists accus- 
tomed to Newtonian mechanics is that while the 
definition of information is subjective (it depends 
on a choice of variables describing the system), 
the processes of destruction and creation of in- 
formation rely on general theories. 

At the risk of being schematic, one may say in- 
formation is created by instabilities (necessarily 
fed by some source of energy) , and it is destroyed 
by relaxation processes (return to equilibrium). 
The resulting stationary patterns are a balance 
between both. A typical laboratory example is 
the Benard Cells. 

In the atmosphere, local hydrodynamical in- 
stabilities result in planetary waves, such as 



the ones responsible for dominant north-westerly 
winds in Canada and south-westerly winds in Eu- 
rope. A pure linear thinker might estimate that 
the hypothetical butterfly that caused the ini- 
tial atmospheric instability is the cause of the 
wave. It is, indeed, chronologically the first of a 
sequence of events that lead to the macroscopic 
pattern. On the other hand, the "non-linear" 
thinker will observe that the macroscopic prop- 
erties of the wave (for example, its spectral char- 
acteristics) do not depend on the position and 
time of the butterfly that triggered the initial in- 
stability. This information has been destroyed by 
the system dynamics. In this view, the "causes" 
of the wave are the conditions that made the ini- 
tial instability possible. 

This lead us to conclude that although there is 
no way to know the climate system fully (it is im- 
possible to know precisely the position and size of 
any molecule of air and ocean and to know all the 
chemical reactions occurring at any given time) 
it is possible to make useful predictions about 
the evolution of some macroscopic variables by 
taking advantage of organised patterns. It is 
therefore sensible to define climate's state using 
variables relevant for what one wants to predict. 
Our goal is to predict the next glacial inception, 
so we will concentrate on variables such as ice 
volume and carbon dioxide concentration. These 
are called order parameters. They describe sys- 
tem's state, but incompletely so, and our goal is 
to establish balance equations for these variables 
(cf. [12 , pp. 36-37). 

Behind the mere desire of predicting the 
next glacial inception, our more general ambi- 
tion is to identify and understand which con- 
straints mostly determine climate evolution at 
the glacial-interglacial time scale. What do we 
need to know to predict the next glacial in- 
ception, and why is this information important 
? This question prompts us to build mod- 
els that take explicitly into account our knowl- 
edge and associated uncertainties, and validate 
these models, that is, to test that the under- 
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lying assumptions are compatible with observa- 
tions (e.g.: [13]). 

We have already seen that instabilities are 
information generators. Macroscopic patterns 
therefore depend on the parameters that con- 
trol the growth of such instabilities. Only in rel- 
atively idealised and simple cases do we know 
these parameters with enough accuracy to cor- 
rectly predict the macroscopic order parameters. 
In most natural cases, instabilities are so numer- 
ous and intricate that the resulting effects cannot 
possibly be predicted without appropriate obser- 
vations. Namely, Saltzman repeatedly insisted 
\14:\ [TB] on the fact that neither current observa- 
tions nor modelling of the present state of the at- 
mosphere can possibly inform us of the ice-sheet 
mass balance with sufficient accuracy to predict 
their evolution at the timescale of several thou- 
sands of years. We need to look at palaeoclimate 
history to get this information. 

The relevant strategy strategy will therefore 
consist in using both first principles and em- 
pirical information to formulate the balance 
equations governing the dynamics of glacial- 
interglacial cycles. Of course, these equations 
must be compatible with our knowledge of atmo- 
sphere and ocean dynamics at the interannual 
time scale, but we accept the fact that we do not 
immediately deduce them from it. The process 
by which model parameters are estimated on the 
basis of observations is called calibration |13j . 

This empirical, or inductive approach is ac- 
ceptable as long as it respects the fundamental 
statements of information theory. In particular, 
any system of time-differential equations that re- 
producibly predicts the evolution of macroscopic 
variables must be dissipative (the volume of ini- 
tial conditions must collapse to an attractor) 
( [TB] and [121, pp. 195 onwards. ) 



2.2 Empirical evidence about the 
Quaternary 

Building a robust theory of glacial-interglacial 
cycles requires a profound knowledge of the Qua- 
ternary. This section is intended to provide a 
glimpse at the vast amount of knowledge that 
scientists have accumulated on that period be- 
fore we propose a mathematical methodology to 
address Ruddiman's hypothesis. 

2.2.1 The natural archives 

By the nineteen- twenties, geomorphologists were 
able to correctly interpret the glacial moraines 
and alluvial terraces as the left-overs of pre- 
vious glacial inceptions. Penk and Briickner 
([T7]. cited by [TH]) recognised four previous 
glacial epochs, named the Giinz, Mindel, Riss 
and Wiirm. The wealth of data on the Quater- 
nary environments that has since been collected 
and analysed by field scientists can be appreci- 
ated from the impressive four-volume encyclo- 
pedia of Quaternary Sciences recently edited by 
Elias [19]. Analysis and interpreting palaeoenvi- 
ronmental data involve a huge variety of scientific 
disciplines, including geochemistry, vulcanology, 
palaeobiology, nuclear physics, stratigraphy, sed- 
imentology, glacial geology and ice-stream mod- 
elling. 

Only a schematic overview of this rich and in- 
tense field of scientific activity could possibly be 
given here. The reader will find most of the rel- 
evant references in the encyclopedia, and only a 
few historical ones are provided here. 

Stable isotopes constitute one important 
class of natural archives. It is indeed known 
since the works of Urey [20] , Buchanan [21] and 
Dansgaard [22 that physical and chemical trans- 
formations involved in the cycles of water and 
carbon fractionate the isotopic composition of 
these elements. To take but a few examples, ice- 
sheet water is depleted in oxygen-18 and deu- 
terium compared to sea water; clouds formed at 
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low temperatures are more depleted in oxygen- 
18 and deuterium than clouds formed at higher 
temperatures; organic matter is depleted in ^^C, 
such that inorganic carbon present in biologically 
active seas and soils is enriched in -"^^C. ^^N is 
another useful stable palaeo-environmental indi- 
cator sensitive to the biological activity of soils. 
The isotopic compositions of water and bio- 
genic carbon are extracted deep-sea sediments, 
ice and air trapped in ice bubbles, palaeosols, 
lake-sediments and stalagmites. One of the first 
continuous deep-sea record of glacial-interglacial 
cycles was published by Cesare Emiliani [23]. 

Radioactive tracers are used to estimate the 
age of the record and rate of ocean water re- 
newal. At the timescale of the Quaternary, useful 
mother-daughter pairs are '^^^Th / 238,234-g (fja.t- 
ing carbonates), and ^"^K / '^'^Ar in potassium- 
bearing minerals. The ratio ^sorjij^ ^ ^^^Pa is a 
useful indicator of ocean circulation rates. 

The chemical composition of fossils is also 
indicative of past environmental conditions. In 
the ocean, cadmium, lithium, barium and zinc 
trapped in the calcite shells of foraminifera indi- 
cate the amount of nutrients at the time of cal- 
cite formation, while the foraminifera content in 
magnesium and strontium are empirically corre- 
lated to water-temperature. 

Glaciologists have also developed ambitious 
programmes to analyse the composition of 
air (oxygen, nitrogen, plus trace gases such as 
methane, carbon-dioxide and nitrogen oxide, ar- 
gon and xenon) trapped in ice accumulating on 
ice sheets, of which the European Project for Ice 
Core in Antarctica is a particularly spectacu- 
lar achievement [24j . It was demonstrated that 
the central plateaus of Antarctica offer a suffi- 
ciently stable environment to reliably preserve 
air's chemical composition over several hundreds 
of thousands of years. The chemical composition 
of water is sensitive to atmospheric circulation 
patterns and sea-ice area. 

Additional sources of informations are ob- 
tained from a variety of marine and continental 



sources. Plant and animal fossils (including pol- 
lens) trapped in lakes, peat-bogs, palaeosols and 
marine sediments provide precious indications on 
the palaeoenvironmental conditions that condi- 
tioned their growth. Their presence (quantified 
by statistical counts) or absence may be inter- 
preted quantitatively to produce palaeoclimatic 
maps [25]. Preservation indicators of ocean cal- 
cite fossils are used to reconstruct the history 
of ocean alkalinity. Palaeosols and wind-blown 
sediments (loess) provide precious indications on 
past aridity at low-latitudes. The loess grain- 
size distribution is also sensitive to atmospheric 
circulation patterns. Geomorphological elements 
remain a premium source of information about 
the configuration of past ice sheets, which is com- 
plemented by datable evidence (typically coral 
fossils) on sea-level. 

2.2.2 The structure of Quaternary cli- 
mate changes 

It is barely straightforward to appreciate which 
fraction of the information available in a cli- 
mate record is relevant to understand climate 
dynamics at the global scale. For example, mi- 
nor shifts in oceanic currents may have sensi- 
ble effects on the local isotopic composition of 
water with however no serious consequence for 
glacial-interglacial cycle dynamics. One strat- 
egy is to collect samples from many areas of the 
world and average them out according to a pro- 
cess called "stacking" . One of the first "stacks" , 
still used today, was published by John Imbrie 
and colleagues [26] in the framework of the Map- 
ping spectral variability in global climate project. 
It is usually referred to as the Specmap stack. 
Here we concentrate on the more recent compila- 
tion provided by Lisiecki and Raymo [27] , called 
LR04. The stack was obtained by superimpos- 
ing 57 records of the oxygen- 18 composition of 
benthic foraminifera shells. Benthic foraminifera 
live in the deep ocean and therefore record the 
isotopic composition of deep water (an indica- 
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Figure 1: The LR04 benthic 6^^0 stack constructed 
by the graphic correlation of 57 globally distributed 
benthic 6^^0 records [27 . Note that the full stack 
goes back in time to -5.2 Myr (1 Myr = 1 million 
years). The signal is the combination of global ice 
volume (low S^^O corresponding to low ice volume) 
and water temperature (low S^^O corresponding to 
high temperature). The Y-axis is reversed as stan- 
dard practice to get "cold" climates down. Data 
downloaded from www.lorraine-lisiecki.com. 

tor of past ice volume). However, there is an 
additional fractionation associated to the calci- 
fication process, which is proportional to water 
temperature. The isotopic composition of cal- 
cite oxygen is reported by a value, named 
giving the relative enrichment of oxygen-18 ver- 
sus oxygen-16 compared to an international stan- 
dard. High 6^^0 indicates either low continental 
ice volume and / or high water-temperature. 

Visual inspection of the LR04 stack (Figure 
[T| nicely evidences the gradual transition from 
the Pliocene — warm and fairly stable — to the 
spectacular oscillations of the late Pleistocene. 
The globally averaged temperature at the early 
Pliocene was about 5° C higher than today ( j28j 
and references therein); that one at the last 
glacial maximum (20,000 years ago) was roughly 
5° C lower. The central research question we 
are busy with is to characterise these oscilla- 
tions, understand their origin and qualify their 
predictability. 



Figure 2: Modulus of the continuous Morlet Trans- 
form of the LR04 stack according to the algorithm 
given in Torrence and Compo [SO] using luq = 5.4 
but with a normalisation c(s) = s/^/At. R routine 
adapted by J.L. Melice and the author from the orig- 
inal code supplied by S. Mallat [91 . The wavelet 
transform evidenes the presence of quasi-periodic sig- 
nals (shades) around periods of 41 kyr (1 kyr = 1,000 
years) and 1000 kyr. 

The Morlet Continuous wavelet transform pro- 
vides us with a first outlook on the backbone of 
these oscillations (Figure [2]) . The LR04 record is 
dominated most of the time by a 40,000-yr sig- 
nal until roughly 900,000 years ago, after which 
the 40,000-yr signal is still present but topped 
by longer cycles. At the very least, this picture 
should convince us that LR04 contains struc- 
tured information susceptible of being modelled 
and possibly predicted. 

How many differential equations will be 
needed? There will be no clear-cut answer 
to that question. Time-series extracted from 
complex systems are sometimes characterised 
by their correlation dimension, which is an 
estimator for the fractal dimension of the 
corresponding attractor |29j . The first estimates 
for the Pleistocene were provided by Nicolis 
and Nicolis [30^ {d = 3.4) and Maasch et al.[3l] 
(4 < d < 6). For this article we calculated 
correlation dimension estimates for the LR04 
stack (d = 1.54) and the HW04 stack g2] 
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{d = 3.56). HW04 is similar to LR04 but it is 
based on different records and dating assump- 
tions. Several authors, including Grassberger 
himself |33l EH ES] have discouraged the use of 
correlation dimension estimates for the "noisy 
and short" time series typical of the Quaternary 
because they are overly sensitive to sampling 
and record length. They are therefore unreliable. 

In response to this problem Ghil and col- 
leagues [351 ES] promoted single-spectrum anal- 
ysis, in which a time series is linearly decom- 
posed into a number of prominent modes (which 
need not be harmonic) , plus a number of small- 
amplitude modes. Assuming that the two groups 
are indeed separated by an amplitude gap, the 
first group provides the low-order backbone of 
the signal dynamics while the second group is 
interpreted as stochastic noise. Single spectrum 
analysis was applied with a certain success to 
various sediment and ice-core records of the few 
last-glacial interglacial cycles [36] and have in 
general confirmed that the backbone of climate 
oscillations may be captured as a linear com- 
bination of a small number of amplitude and 
/ or frequency-modulated oscillations. Single- 
spectrum analysis of the last million years of 
LR04 (Figure |3]) confirms this statement. 

2.2.3 The Achille heel 

Now time has come to mention a particularly dif- 
ficult and intricate issue: dating uncertainty in 
palaeoclimate records. No palaeoclimate record 
is dated with absolute confidence. Marine sedi- 
ments are coarsely dated by identification of a 
number of reversals of Earth's magnetic field, 
which have been previously dated in rocks by 
radiometric means ([37] and references therein). 
Magnetic reversals are pretty rare (four of them 
over the last 3 million years) and their age is 
known with a precision no better than 5,000 
years. Local sedimentation rates may vary con- 
siderably between these time markers such that 
any individual event observed in any core taken 



in isolation is hard to date. Irregularities in the 
sedimentation rate blur and destroy information 
that might otherwise be evidenced by spectral 
analysis. 

One strategy to contend this issue is to assume 
synchrony between oscillation patterns identi- 
fied in different cores. Statistical tests may 
then be developed on the basis that dating er- 
rors of the different cores are independent. For 
example, Huybers (2007) [38] considered the 
null-hypothesis that glacial-interglacial transi- 
tions (they are called terminations in the jar- 
gon of palaeoclimatologists) are independent on 
the phase of Earth's obliquity. While this null- 
hypothesis could not be rejected on the basis of 
a single record, the combination of 14 cores al- 
lowed him to reject it with 99% confidence, prov- 
ing once more the effect of the astronomical forc- 
ing on climate. First tests of this kind were car- 
ried out by Hays et al. in a seminal paper [39j. 
Note that in many cases the oscillation patterns 
recognised in different cores are so similar that it 
is hard to dispute the idea of somehow "matching 
them", but it is remarkable that rigorous statis- 
tical tests assessing the significance of a correla- 
tion between two ill-dated palaeoclimate records 
are only being developed (Haam and Huybers, 
manuscript in preparation). 

Another strategy is known as orbital tuning. 
The method consists in squeezing or stretching 
the time-axis of the record to match the evolu- 
tion of one or a combination of orbital elements, 
possibly pre- filtered by a climate model |39| |26] . 
The method undeniably engendered important 
and useful results (e.g. [ID]), but the astute 
reader has already perceived its potential per- 
versity: orbital tuning injects a presumed link 
between orbital forcing and the record. Experi- 
enced investigators recognise that orbital tuning 
has somehow contaminated most of most of the 
dated palaeoclimate records available in public 
databases. This has increased the risk of tauto- 
logical reasoning. 

For example, compare the two SSA analyses 
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Figure 3: Single Spectrum Analysis (SSA) of the LR04 and HW04 benthic stacks. Displayed are the eigenvalues 
of the lagged-covariance matrix of rank AI = 100 as given by [92^, equation (6). The records were cubic-spline 
interpolated {At — Ikyr) and only the most recent 900 kyr were kept. The SSA decomposition of LR04 is very 
typical: it evidences three oscillators (recognisable as pairs of eigenvectors), then about four modes that are 
generally interpreted as harmonics of the dominant ones, and finally a number of modes typically interpreted 
as stochastic background. The HW04 stack contrasts with LR04 because the dominant modes are not so easily 
evidenced. HW04 uses less benthic records than LR04, but it also relies on more conservative dating assumptions 
and this probably resulted in blurring the quasi-periodic components of the signal. HW04 data were obtained 
from www . people . f as . harvard . edu/ phuybers/. 



shown in figure [3| As we mentioned, LR04 and 
HW04 are two stacks of the Pleistocene but LR04 
contains more information. It is made of more 
records (57 instead of 21 in HW04) and it is as- 
tronomicahy tuned. We can see from the SSA 
analysis that LR04 presents more quasi-periodic 
structures than HW04 (recall that quasi-periodic 
modes are identified as pairs of eigenvalues with 
quasi the same amplitude). Why is this the case? 
Is this because age errors in IIW04 blurred the 
interesting information, or is it because this in- 
formation has been artificially injected in LR04 
by the tuning process? There is probably a bit 
of both (but note that IIW04 displays a similar 
wavelet structure as LR04). 

Leads and lags between CO2 and ice volume 
is another difficult problem where risks posed by 
hidden dating assumptions and circular reason- 
ing lie at every corner. Here is one typical illus- 
tration: Saltzman and Verbitsky showed at sev- 
eral occasions (e.g.: [^P) a phase diagram show- 
ing the SPECMAP 5^^0 stack versus the first 



full ice-core records of CO2 available at that time 
[i2l Ii3] . It is reproduced here (Figure |4] left). 
The phase diagram clearly suggests that CO2 
leads ice volume at the 100-kyr time scale. How- 
ever, a detailed inspection of the original publi- 
cations reveals that the SPECMAP record was 
astronomically tuned, and that the Vostok time- 
scale uses a conventional date of isotopic stage 
5.4 of 110 kyr BP ... by reference to SPECMAP 
|43] ! The hysteresis is therefore partly condi- 
tioned by arbitrary choices. On Figure |4] we fur- 
ther illustrate the fact that the shape of the hys- 
teresis depends on the stack record itself. The 
situation today is that there is no clear consen- 
sus about the phase relationship between ice vol- 
ume and CO2 at the glacial interglacial time scale 
(compare [SI 051 SB]). According to the quite 
careful analysis of [l5], CO2 leads ice volume at 
the precession (20 kyr) period, but CO2 and ice 
volume are roughly synchronous at the obliquity 
(40 kyr) period. Current evidence about the lat- 
est termination is that decrease in ice volume 
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and the rise in CO2 were grossly simultaneous 
and began around 19,000 years ago |44| ITT] 

2.3 Getting physical laws into the 
model 

So far we learned that palaeoclimate oscillations 
are structured and that it is not unreasonable 
to attempt modelling them with a reduced order 
model forced by the astronomical variations of 
Earth's orbit. What is the nature of the physi- 
cal principles to be embedded in such a model, 
and how can they be formalised? The history 
of Quaternary modelling is particularly enlight- 
ening in this respect (the reader will find in [18j 
an extensively documented review of Quaternary 
climates modelling up to the mid-eighties). Af- 
ter Joseph Adhemar (1797 - 1862) [48 suggested 
that the cause of glaciations is the precession of 
the equinoxes, Joseph John Murphy and James 
CroU (1821 - 1890) argued about how preces- 
sion may affect climate. Murphy maintained 
that cold summers (occurring when summer is 
at aphelion) favour glaciation while Croll 
considered that cold winters are critical |50] . 

Croll's book demonstrates a phenomenal en- 
cyclopaedic knowledge. His judgements are 
at places particularly far-sighted, but they are 
barely substantiated by the mathematical analy- 
sis Fourier was so much insistent about. The na- 
ture of his arguments are essentially phenomeno- 
logical, if not at places frankly rhetorical. 

Milutin Milankovitch (1879 - 1958) is then 
generally quoted as the one having most deci- 
sively crossed the step towards mathematical cli- 
matology. In a highly mathematical book that 
crowns a series of articles written between 1920 
and 1941 [51], Milankovitch extends Fourier's 
work to estimate the zonal distribution of Earth's 
temperature from incoming solar radiation. He 
also computes the effects of changes in preces- 
sion, eccentricity and obliquity on incoming solar 
radiation at different latitudes to conclude, based 
on geological evidence, that summer insolation is 



indeed driving glacial-interglacial cycles. 

Mathematical analysis is the process that al- 
lows Milankovitch to deduce the consequences of 
certain fundamental principles, such as the laws 
of Beer, Kirchhoff and Stefan, on global quanti- 
ties such as Earth's temperature. On the other 
hand, Milankovitch uses empirical macroscopic 
information, such as the present distribution of 
the snow-line altitude versus latitude, to esti- 
mate the effects of temperature changes on the 
snow cover. In today's language, one may say 
that Milankovitch had accepted that some infor- 
mations cannot be immediately inferred from mi- 
croscopic principles because they depend on the 
way the system as a whole has been dealing with 
its numerous and intricate constraints (Earth's 
rotation, topography, air composition etc.). 

The marine- record study published by Hays, 
Imbrie and Shackleton [39^ is often cited as the 
most indisputable proof of Milankovitch 's theory. 
Hays et al. identified three peaks in the spectral 
estimate of climate variations that precisely cor- 
respond to the periods of obliquity (40 kyr) and 
precession (23 kyr and 19 kyr) calculated analyt- 
ically by Andre Berger|^ 

However, sensu stricto, Milankovitch's theory 
of ice ages was invalidated by evidence — al- 
ready available in an article by Broecker and van 
Donck [56] — that the glacial cycle is 100,000 
years long, ice build up taking about 80,000 years 
and termination about 20,000 years |56| l39] . 
Neither the 100,000-year duration of ice ages, 
nor their saw-tooth-shape were predicted by Mi- 
lankovitch. The bit Milankovitch's theory is 
missing is the dynamical aspect of climate's re- 
sponse. Glaciologist Weertman [57j consequently 
addressed the evolution of ice sheet size and vol- 
ume by means of an ordinary differential equa- 
tion, thereby opening the door to the use of dy- 
namical system theory for understanding Qua- 

^the supporting papers by Berger would only appear 
in the two following years [521 1531 154j ; Hays et al. based 
themselves on a numerical spectrum estimate of the or- 
bital time-series provided by Vernekar [55] 
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Figure 4: The concentration in CO2 measured in the Vostok ice core record fSD] over the last glacial-interglacial 
cycle is plotted versus two proxies of continental ice volume: (left) : The planctonic S^^O stack by Imbrie et al. 
(1984) and (right) : The benthic S^^O stack by Lisiecki and Raymo (2004). Numbers are dates, expressed in 
kyr BP (before present). While the Imbrie stack suggests an hysteresis behaviour with CO2 leading ice- volume 
variations, the picture based on LR04 is not so obvious. 



ternary oscillations. 

In the meantime, general circulation models 
of the atmosphere and oceans running on su- 
percomputers became widely available (cf. |58j 
for a review), and used for palaeoclimate pur- 
poses [591 EQI E] • The interest of these models 
is that they provide a consistent picture of the 
planetary dynamics of the atmosphere and the 
oceans. Just as Milankovitch applied Beer and 
Kirschoff 's laws to infer Earth's temperature dis- 
tribution, general circulation models allow us to 
deduce certain aspects of the global circulation 
from our knowledge of balance equations in each 
grid cell. However, these balance equations are 
uncertain and quantifying the consequences of 
these uncertainties at the Earth global scale is a 
very deep problem that only begins to be system- 
atically addressed ^62^ . While general circulation 
models are undeniably useful to constrain the im- 
mediate atmospheric response to changes in or- 
bital parameters, they are far too uncertain to 
reliably estimate glacial accumulation rates with 



enough accuracy to predict the evolution of ice 
sheets over tens of thousands of years |14j . 

In the following sections we will concentrate on 
a 3-dimensional climate dynamical model writ- 
ten by Saltzman. This choice was guided by 
the ease of implementation as well as the im- 
pressive amount of supporting documentation 
|15] . However, they were numerous alterna- 
tives to this choice. The reader is referred to 
the article by Imbrie et al. [B^ and pp. 264- 
265 of Saltzman's book [15^ for an outlook with 
numerous references organised around the dy- 
namical concepts proposed to explain glacial- 
interglacial cycles (linear models, with or with- 
out self-sustained oscillations, stochastic reso- 
nance, model with large numbers of degrees of 
freedom) . 

The series of models published by Ghil and col- 
leagues [Ml ESI [66] are among the ones having the 
richest dynamics. They present self-sustained os- 
cillations with a relatively short period (6,000 
years). The effects of the orbital forcing are 
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taken into account by means of a multiplicative 
coefficient in the ice mass balance equation. This 
causes non-linear resonance between the model 
dynamics and the orbital forcing. The result- 
gin spectral response presents a rich background 
with multiple harmonics and band-limited chaos. 
More recently, Gildor and Tziperman [67J pro- 
posed a model where sea-ice cover plays a cen- 
tral role. In this model, termination occurs when 
extensive sea-ice cover reduces ice accumulation 
over ice sheets. Like Saltzman's, this model 
presents 100-kyr self-sustained oscillations that 
can be phase-locked to the orbital forcing. 

Field scientists with life-long field experience 
have also proposed models usually qualified as 
"conceptual", in the sense that they are formu- 
lated as a worded causal chain inferred from a 
detailed inspection of palaeoclimate data with- 
out the support of differential equations. Good 
examples are |l5l [631 EH EH]- In the two lat- 
ter references, Ruddiman proposes a direct effect 
of precession on CO2 concentration and tropi- 
cal and southern-hemisphere sea-surface temper- 
atures, while obliquity mainly affects the hydro- 
logical cycle and the mass-balance of northern 
ice sheets. 

2.4 The Saltzman model (SM91) 

As a student of Edward Lorenz, Barry Saltz- 
man ( - 2002) contributed to the formulation and 
study of the famous Lorenz63 dynamical system 
[70] traditionally quoted as the archetype of low- 
order chaotic system]^ Saltzman was therefore in 
an excellent position to appreciate the explana- 
tory power of dynamical system theory. Between 
1982 and 2002 he and his students published a 
small dozen of dynamical systems deemed to cap- 
ture and explain the dynamics of Quaternary os- 
ciUations [HI O EH Ei [E] . In the present 

^The acknowledgments of the Lorenz (1963) paper 
reads: "The writer is indebted to Dr. Barry Saltzman 
for bringing to his attention the existence of nonperiodic 
solutions of the convection equations. 



article we choose to analyse the "palaeoclimate 
dynamical model" published by Saltzman and 
Maasch (1991) [73 . We will refer to this model 
as SM91. 

Saltzman estimated that the essence of Qua- 
ternary dynamics should be captured by a three- 
degree-of-freedom dynamical system, possibly 
forced by the variations in insolation caused by 
the changes in orbital elements [13]. The evolu- 
tion of climate at these time scales is therefore 
represented by a trajectory in a 3-dimensional 
manifold, which Saltzman called the "central 
manifold". The three variables are ice vol- 
ume (/), atmospheric CO2 concentration (/i) and 
deep-ocean temperature {6). It is important to 
realise that Saltzman did not ignore the existence 
of climate dynamics at shorter and longer time 
scales than those that characterise the central 
manifold, but he formulated the hypothesis that 
these modes of variability may be represented by 
distinct dynamical systems. In this approach, 
the fast relaxing modes of the complex climate 
system are in thermal equilibrium with its slow 
and unstable dynamical modes. This assump- 
tion is called the "slaving principle" and it was 
introduced by Haken [75] . 

The justification of time-scale decoupling is a 
very delicate one and it deserves a small digres- 
sion. In some dynamical systems, even small 
scale features may truly be informative to pre- 
dict large-scale dynamics. This phenomenon, 
called "long-range interaction", happens in the 
Lorenz63 model [76]. The consequence is that 
one might effectively ignore crucial information 
by averaging the fast modes and simply assume 
that they are in thermal equilibrium. To justify 
his model, Saltzman used the fact that there is 
a "spectral gap" , that is a range of periods with 
relatively little variability, between weather (up 
to decadal time-scales) and climate (above one 
thousand years). This gap indicates the pres- 
ence of dissipative processes that act as a barrier 
between the fast and the slow dynamics. It is 
therefore reasonable to apply the slaving princi- 
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pie. In relation to this, Huybers and Curry re- 
cently published a composite spectral estimate of 
temperature variations ranging from sub-daily to 
Milankovitch time scales [T^. No gap is evident, 
but Huybers and Curry identify a change in the 
power-law exponent of the spectral background: 
Signal energy decays faster with frequency at the 
above the century time scale than below. They 
interpret this as an indice that the effective dis- 
sipation time scale is effectively larger above the 
century that below, and, therefore, that the dy- 
namics of slow and fast climatic oscillations are 
at least partly decoupled. 

We now enunciate the three differential equa- 
tions of SM91. 

The ice-mass balance is the result of the con- 
tribution of four terms: a drift, a term inversely 
proportional to the deviation of the mean global 
temperature compared to today (r), a relax- 
ation term, and a stochastic forcing representing 
"all aperiodic phenomena not adequately param- 
eterised by the first three terms" : 

^ = ifi- ^2f- ^3! + Wiit). (1) 
at 

According to the slaving principle, r is in ther- 
mal equilibrium with the slow variables {/, fi, 9} 
and its mean may therefore be estimated as a 
function of the latter: 

f = friI)+f^,ifi)+fe{9)+fR{R), (2) 

where fx{x) is the contribution variation of x 
compared to a reference state, to f keeping the 
other slow variables or forcing constant. R des- 
ignates the astronomical forcing (Saltzman used 
incoming insolation at 65° N at summer solstice). 
The different terms r(.) were replaced by linear 
approximations, the coefficients of which were es- 
timated from general circulation model experi- 
ments. 

The CO2 equation includes the effects of ocean 
outgassing as temperature increases, a forcing 
term representing the net balance of CO2 in- 
jected in the atmosphere minus that eliminated 



by silicate weathering, a non-linear dissipative 
term and a stochastic forcing: 

with Kfj_ = f3i- [32^1 + iSsn'^. 

The dissipative term {K^fj,) is a so-called Landau 
form and its injection into the CO2 equation is 
intentional to cause instability in the system. In 
an earlier paper (e.g. [ZH])) Saltzman and Maash 
attempted to justify similar forms for the CO2 
equation on a reductionist basis: each term of 
the equation was identified to specific, quantifi- 
able mechanisms like effect of sea-ice cover on 
the exchanges of CO2 between the ocean and 
the atmosphere or that of the ocean circulation 
on nutrient pumping. It is telling that Saltzman 
and Maash gradually dropped and add terms to 
this equation (compare [7S1 [721 US] ) to arrive at 
the above formulation by which they essentially 
posits a carbon cycle instability without explic- 
itly caring about which mechanisms that cause 
it. 

The deep-ocean temperature simply assumes 
a negative dependency on ice volume with a dis- 
sipative relaxation term: 

d6 

-^=li-l2l- + We (4) 
at 

The carbon cycle forcing term Ffj_ is assumed 
to vary slowly at the scale of Quaternary os- 
cUlations. It may therefore be considered to 
be constant and its value is estimated assuming 
that that the associated equilibrium is achieved 
for a CO2 concentration of 253 ppmv. We 
shall note {/o,/Uo,^o} the point of the central 
manifold corresponding to that equilibrium, and 
{I',fi',6'} the departure from it. Further con- 
straints are imposed by semi-empirical knowl- 
edge on the relaxation times of ice sheet mass 
balance (10,000 years) and deep-ocean tempera- 
ture (4,000 years). 
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Saltzman and Maasch explored the different 
solution regimes of this system \79\ [73] and they 
observed that climate trajectories converged to 
a limit cycle characterised by saw-tooth-shaped 
oscillations for a realistic range of parameters. 
When the model is further forced by the astro- 
nomical forcing, the uncertainty left on the em- 
pirical parameters of equations ([T] - 14]) provides 
the freedom to obtain very convincing solutions 
for the variations in ice volume and CO2 dur- 
ing the late Quaternary. Figure [5] reproduces 
the original solution [73J, using the parameters 
published at the time. As in the original pub- 
lication, the solution is compared with Imbrie's 
(5^^0-stack [26 interpreted as a proxy for ice vol- 
ume, and CO2 record extracted from the Vostok 
and EPICA(Antarctica) ice cores [801 l8T]. 



Dynamical model of SM91 



Limit-cycle solutions in SM91 (Figure 2.4) owe 



their existence to cubic terms in the CO2 equa- 
tion. In fact, all parameters being constant, a 
limit-cycle occurs only for certain carefully cho- 
sen values of fiQ, which led Saltzman to conclude 
that the cause of glacial-interglacial oscillations 
is not the astronomical forcing (a linear view of 
causality) but rather the gradual draw-down of 
/io at the tectonic time scale that permitted the 
transition between a stable regime to a limit- 
cycle via a Hopf bifurcation. According to this 
approach, astronomical forcing controls in part 
the timing of terminations by a phase-locking 
process, but terminations essentially occur be- 
cause negative feedbacks associated to the car- 
bon cycle become dominant at low CO2 concen- 
tration and eject the system back towards the 
opposite region of its phase space. 
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Figure 5: Response of the palaeoclimate model of 
Saltzman and Maasch (1991) [73]. Shown are the 
insolation forcing, taken as the summer solstice in- 
coming solar radiation at 65° N after [54]; the ice 
volume anomaly (full), overlain with the SPECMAP 
planctonic S^^Oc stack ^ (dashed), the CO2 atmo- 
spheric concentration, overlain with the Antarctic ice 
core data from Vostok and EPICA [101 [H], and fi- 
nally deep-ocean temperature. Note that /' and 9' 
are anomalies to the tectonic average. A similar fig- 
ure was shown in the original article by Saltzman and 
Maasch 



3 The Bayesian inference pro- 
cess 

Approaches founded on low-order dynamical sys- 
tems are regularly suspected of being tautolog- 
ical: what can you learn from a model if you 
tuned it to match observations? There is no 



doubt that the empirical content of any model 
— i.e., its capacity of being in conflict with ob- 
servations — has to be assessed with the utmost 
care. Several authors have in particular insisted 
on the difficulty of finding discriminating tests 
for models with similar dynamical characteris- 
tics but built on different interpretations of the 
climate system's functioning ^82^ i83j . It is there- 
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(a) SM91 model without astronomical forcing 



(b) SM91 model with astronomical forcing 
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Figure 6: Phase-space diagrams 
of trajectories simulated with the 
SM91 model, using standard pa- 
rameters. The model exhibits a 
limit cycle in absence of external 
forcing, with a trajectory that re- 
sembles those obtained with data 
(Figure [4]). The astronomical 
forcing adds a number of degrees 
of freedom that complicates the 
appearance of the phase diagram. 



fore challenging but important to identify and 
design powerful tests for such models. 

Nevertheless, it has to be appreciated that 
the risk of tautology is also present in the most 
sophisticated general circulation models of the 
atmosphere and ocean. Once assembled, these 
models are "tuned" to capture major and global 
characteristics of climate such as the overturning 
cell or the global mean temperature (e.g.: [Hi]). 
This "tuning" is an effective way of incorporating 
macroscopic information in the model, and this 
information can no longer said to be "predicted" . 

Statistical decision theory allows us to address, 
at least partly, these difficult problems. We will 
concentrate on one branch of it: Bayesian infer- 
ence. The paradigm of Bayesian inference finds 
its roots in early works by Bayes, Laplace and 
Bernouilli who were looking for ways of augment- 
ing their knowledge of certain quantities such 
as initial conditions or parameters, by means of 
observations [16 . Rougier (2007) [13] explains 
how Bayesian inference methods may be applied 
to the problem of climate prediction. His con- 
clusions are summarised hereafter, but adapted 
were relevant to the problems posed by palaeocli- 
mate time-series analysis. Compared to Rougier 
(2007), we more explicitly consider here the fact 
that climate is a dynamical body, whose evo- 
lution has to be predicted by means of time- 
differential equations. 

Before embarking on the mathematical de- 
tails, it is useful to recall two aspects inherent 



to complex system modelling introduced in sec- 
tion 2.1 The first one is that by focussing on 
certain modes of climate variability we ignore a 
large body of information, such as its synoptic 
variability and, for example, the occurrence of 
a particular volcanic eruption at any particular 
moment. This ignorance causes prediction er- 
rors that we have to parameterise, typically as a 
stochastic forcing or error (we will here neglect 
the epistemological distinction between stochas- 
tic forcing and error). Model validation con- 
sists in verifying that the model assumptions are 
compatible with observations. A crucial but of- 
ten forgotten point is that the validation tests 
depends on the judgements we will have made 
about the probability distribution of the model 
error: if we considered that the model error could 
take any value, the model would always be com- 
patible with observations, but it would also be 
useless. 

The second aspect of complex system mod- 
elling is that we accept to consider information 
that it is not immediately deduced from our 
knowledge of microscopical interactions. In the 
case we are busy with, these extra statements 
take the form of conjectures about the math- 
ematical expressions of carbon, ocean and ice- 
sheet feedbacks, which are calibrated by reference 
to observations. 

Our purpose is to formalise as rigorously as 
possible the validation and calibration processes. 
To this end, let us denote y{t) a vector describing 
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the state of climate at a given time t. We further 
notate symbolically y the climate evolution over 
a given time interval not necessarily restricted 
to the observable past. It is useful to distinguish 
notationally the variable Y, which may a priori 
take any value in a given space, from its realisa- 
tion y. The exact value of y is never known be- 
cause any measurement or prediction is affected 
by errors, but the fact of positing the existence 
and attaching a meaning to y enables us to struc- 
ture and justify our judgements. 

Palaeoclimatologists attempt to retrieve in- 
formation on y by taking measurements in a 
palaeoenvironmental record. Let z he a series 
of such observations like, for example, delta- 
Deuterium of ice in an Antarctic ice core sampled 
at certain depths. They estimate that z is con- 
ditionally dependent on y, which one may write 
as : 

y z (5) 



This means that their expectation on z depends 
on their knowledge of y. This expectation can 
be quantified by means of a probability density 
function for Z, thought of as a function of z: 

P{Z = z\Y = y,p) (6) 

Building an expression for ^ requires to for- 
mulate a number of assumptions forming a cli- 
mate proxy model that we have symbolically de- 
noted p. In practice, it may be preferable to de- 
compose this model into a chronological chain of 
nested processes, each bearing uncertainties: ef- 
fect of climate on the hydrological cycle, isotopic 
fractionation, accumulation of ice, preservation 
of the signal in the core, drilling and actual mea- 
surement. The more there are uncertainties, the 
wider P(Z = z\Y = y,p) will be. 

Bayesian inversion then indicates us how z is 
informative on y : 



P{Y = y\Z = z,p) 



P{Z = z\Y = y)P{Y = y) 



P{Z = z\p) 



(7) 



This equation bears important lessons. First, up- 
dating an estimate of y on the basis of obser- 
vations requires to have some prior judgement 
expressed in the form P(Y = y). This im- 
portant question will be kept aside a moment. 
Second, the denominator at the right-hand-side 
is independent on y. It represents a marginal 
likelihood, which may be thought of as a point- 
estimate of a predictive distribution of Z given 
our prior judgement on y along with the assump- 
tions contained in p. In practice it is evaluated 
as: 



P{Z = z\p) = j P{Z = z\Y = y,p)P{Y = y) dy 

(8) 

The validation of p consists in determining if 
P{Z = z\p) lies in the tails of its distribution. 
The presence of an observation in the tails of its 
predictive distribution means that it was little 
likely to occur according to the theory expressed 
in p. Such an outcome will incline us to con- 
fidently reject the theory in the same way that 
one rejects a null-hypothesis in classical statistics 
tests. This is easily diagnosed in the case where 
z is a scalar, in which case it may be checked if 
the marginal probability P{Z < z\p) is not too 
close to zero or one. 

P{Z < z\p) = j P{Z < z\Y = y,p)P{Y = y) dy 

(9) 

In practice, z is often highly dimensional and its 
predictive distribution may be particularly intri- 
cate, especially in chaotic dynamical systems. 

At present, it is useful to split y into its "past" 
(yp) and "future" (y/) components. If the past 
is known, the record content is obviously inde- 
pendent on the future, i.e. : 

P{Z = z\Yp = yp,Yf = yj) = P{Z = z\Yp = y^). 

(10) 

([T]) and (10) tell us that in absence of any 
additional assumption, past observations 
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are not informative on the future. Predict- 
ing climate requires to assume a certain dynami- 
cal structure to climate evolution to link yj to i/h- 
This is the role of the climate model. It is in prin- 
ciple always possible to formulate this model in 
terms of first-order differential stochastic equa- 
tions if the climate state y{t) is suitably defined. 
Climate time-series are in this case Markovian: 
Given climate at any time to, the probability 
density function of climate at time ti may be 
estimated and written: 



P{yihMto),c,A = a) 



where we distinguish the ensemble of model 
equations (symbolically denoted c) from their pa- 
rameters, gathered into a single vector variable 
denoted A. More generally, the model allows us 
to estimate the probability density of any climate 
time-series, which we shall write: 



y{to) 



(12) 



The model makes it thus possible to build a 
predictive distribution function for y given any 
prior estimate of the possible values of a and 

y(o). 

The climate and climate proxy models may 
then be combined to form a Bayesian network 



y{to) 



■Vh 



(13) 




Solving the network means to find the joint 
distribution of a, y(to)) Vh-, Up and z compatible 
with all the constraints expressed in p and c (e.g. 
: [85], pp. 167 and onwards). Keeping in mind 
that the arrows may be "inverted" by application 



of Bayes' theorem, it appears that there are two 
routes by which z constrains yp : via yh (that 
is, constraining the initial conditions to be input 
to the model forecast of the future), and more 
indirectly via a. In the latter route, all obser- 
vations concur to constrain a distribution of the 
model parameters that is compatible with both 
the model structure and the data. 



Two more remarks. First, (13) shows that the 



climate model has solved the problem of finding 
a prior to y (it is provided by the model), but 
this is at the price of having to find a prior for 
parameters a. It may happen that one parameter 
has no clearly identified physical meaning (like /54 
in eq. ([3|) and we would like to express our total 
ignorance about it, except for the fact that it is 
positive. It happens that there is no definitive 
solution to the problem of formulating a totally 
ignorant prior. However, if the observations are 
very informative, the posterior distribution of a 
is expected to depend little on its prior. 

The second remark is about the marginal like- 
lihood, that is, our assessment of the plausible 
character of observations z given the structural 
assumptions in models p and c along with the 
prior on a. It is crucial to be clear about what 
is being tested. For example, one may be con- 
tent to assess the position of z thought of as a 
n-dimensional vector (n is the number of obser- 
vations) in the manifold of likely Z values given 
the prior on a. This test takes for granted that 
the stochastic error is effectively white-noise dis- 
tributed. This being said, it may be useful to 
effectively test the white-noise character of the 
model error, typically example by estimating the 
likelihood of the lagged-correlation coefficients 
of the stochastic error. Lagged-correlation co- 
efficients significantly different than zero almost 
surely indicate that there is information in the 
stochastic error terms. This would prove that 
the model is incomplete in the sense that its pre- 
dictive ability can almost surely be improved. 
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An application of the particle 
filter 



The important difference with (13) is that the 



^Network ^ is an example of combined param- 
eter and time-varying state estimation problem. 
This kind of problem is highly intractable, but 
statisticians have been looking at ways of find- 
ing approximate solutions based on Monte-Carlo 
simulations. Here we use an implementation of 
the particle filter developed by Liu and West |86] . 
This is a filter, that is a sequential assimilation 
method: observations are used to refine parame- 
ter distribution estimates as the time-integration 
of the model progresses. The reader is referred 
to the original publication for a fuller discussion 
of the method and we will briefly summarise here 
the sequential algorithm. 



First we reformulate ( 13 ) into a more tractable 
problem: 



(14) 




*This section outlines work in progress carried out by 
the author in close collaboration with Jonathan Rougier, 
Department of Statistics at the University of Bristol, UK. 



observations are bound to individual state vec- 
tors. This implies that their dating is certain 
(they can unambiguously be associated to a cli- 
mate state at a given time) and that there is no 
diffusion of the signal within the record. 

The climate model (c) is SM91 ([l]-!!]), the 
equations of which are summarised hereafter: 



dl^ 
dt 

d4j!_ 
~dt 

~dt 



-ai[kf,n' + keO' + knR'it)] - Kjl' + W/ 

(cl) 

hill' - 62/x'2 + 63/^'^ - bgO + (c2) 



(c3) 



The coefficients Oj, 6j, Cj and Ki are functions 
of the 7i determined using the condition 

that the equations for {I' , fi' ,9'} present a fixed- 
point at (i.e., {/o,/xoi^o} is a long-term, "tec- 
tonic" equilibrium). Coefficients kx appear in 
the process of linearising the short-term response 
and can in principle be estimated with general 
circulation models. The reader is referred to the 
original publications for fuller details. 

The climate proxy model {p) is very simple. 
We will use the Specmap stack of planctonic 
foraminifera to constrain ice volume [26 , and 
the Vostok (Antarctica) ice core by Petit et al. 
(1999) [80_ for CO2. 



CO2 



0.71 

45 IQis m3 
fJ- + Wco2 



(Pl) 
(p2) 



Equation (pl) uses the fact that the Imbrie 
et al. record is expressed in standard deviation 
units with zero mean, along with the constraint 
that a total ice melt of 45 10^^ m^ is recorded 
as a drop of 0.71 (unitless) in Imbrie et al. We 
therefore neglect the influence of ocean tempera- 
ture on the record, while this issue is contentious. 
Errors are parameterised by means of additive 
stochastic Gaussian white noise with standard 
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deviations of 0.2 (pi) and 20 ppm (p2), respec- 
tively. 

The above approximations (neglecting dating 
uncertainty, in-core diffusion and unduly simple 
isotope model) will no longer be tenable as this 
research project develops but they are suitable 
for a first application of the particle filter algo- 
rithm. Consequently, results should be consid- 
ered with the necessary caution. 

We now review the particle-filter algorithm. A 
particle is essentially a realisation of the state 
vector (say : 2/(to)) associated to a realisation of 
the parameters {A = {ln(aj, 6j, q)}) and a weight 
{w). Ten thousand (n) particles are initialised 
by sampling the prior of y{tO) and A. Prior pa- 
rameter distribution are log-normal around the 
values given in SM91 (Figure [7]). Only the Oj, bi, 
ci and kg are considered to be uncertain, while 
the dissipative exchange coefficients Kj and Kg 
as well as the climate sensitivities and kn are 
assumed to be known (Table [T]). 

All weights are initialised to 1. The filter then 
consists of an iterative six-step process. Say we 
are at time t. 

1. Propagation, that is, time- integration of 
all particles until the time (t+l) correspond- 
ing to the next available data (either CO2 or 

2. Shrinkage. Particles are now dispersed in a 
region of the {Y, A}. This region in shrunk, 
that is, the particles are made closer to each 
other by a factor a. 

3. Weight estimate Particle weights are mul- 
tiplied by the likelihoods P{Z = z\Y = yj), 
where yj is the state of particle j, and z is 
the encountered data. 

4. Importance resamphng based on pos- 
terior estimate. After step 2, some parti- 
cles may be given a large weight while others 
only a small one. Particles are therefore re- 
sampled in such a way that they all get a 



a, (m^/K/yr) 



b, (1/yr) 



b, (1/ppm/yr) 



b, (l/ppm'/yr) 





5e-07 3e-06 



bj (ppm/K/yr) 



b, (K/mVyr) 





0.001 0.004 0.015 



kg (unitless) 



2e-24 8e-24 




Figure 7: Prior (dashed) and posterior (full) density 
estimates of the parameters allowed to vary in SM91. 
The filter has been successful in narrowing down the 
distributions. 

similar weight. This implies that some par- 
ticles are duplicated while others are killed. 
Particles are now distributed along k < n 
kernels. 

5. Resamphng of kernels Each kernel is bro- 
ken apart into particles with parameters 
scattered with variance /i^. 

6. Weight update Particles weights are up- 
dated according to their likelihood. 

Shrinkage and kernel sampling are artefacts in- 
troduced to avoid filter degeneracy. Liu and 
West [86] note that the estimator is unbiased for 
a = (3(5 - l)/2d and h"^ = 1 - o? . The param- 
eter b is called a discount factor. It must lie in 
]0, 1] and typically around 0.95 - 0.99. Here we 
choose 5 = 0.95. We found that the parameter 
disturbance due to the filter dominates any rea- 
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4 AN APPLICATION OF THE PARTICLE FILTER 



IJ 240 
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Figure 8: Filtered state estimates with the SM91 
model constrained by the SPECMAP data (squres), 
and the Antarctic ice core data (pluses). The state 
estimates are represented by shades, dark and light 
gray representing the [25**^; 75"^] and [5**^; 95"^] quan- 
tiles of the particle weighted distributions, respec- 
tively. The lower graph represents, for each data, 
the model predictive probability that the data would 
have been lower than it actually was, given the pre- 
vious parameter and state estimates. The repetition 
of probabilities below 0.05 or above 0.95 tend to in- 
validate the model. 

sonable amount of stochastic error that could be 
parameterised via W. Therefore, we decided not 
to account for the model stochastic error noise 
to gain computing efficiency. 

Figure [8] summarises the essential features of 
the particle filter run. It represents, for each 
prognostic variable, the evolution of the state es- 
timate (shaded) along with the data. The dark 
and light shades represent the central 50 and 90 
% percentiles of the weighted particle distribu- 
tion. The filter algorithm updates the parameter 
estimates as it meets the data (The posterior pa- 
rameter distributions are compared to the prior 



on Figure [7]), which explains why the state esti- 
mates become narrowed as time progresses. The 
dots and pluses are the observation estimates of 
ice volume and CO2. The fourth panel is a first 
step towards model validation. It display, for 
each observation, the model predictive probabil- 
ity that this observation was smaller or equal 
than its value, exactly in the spirit of equation 
([9]). Values too close to zero or one cast doubt 
on the model. 

It was unexpected that the fit of the state esti- 
mates of the ice volume on SPECMAP would be 
so poor. In fact, the model systematically over- 
estimates ice volume during interglacials and this 
occurs as soon as CO2 observations are taken into 
account. Strictly speaking, the model is invali- 
dated. Where does the problem lie ? 

The most obvious possibility is that we have 
incompletely modelled the SPECMAP stack. In- 
deed, we know that water temperature con- 
tributes to the S^^Oc signal but this contribu- 
tion is missing in the model ( |i6 l 187 1 188 | |89] . the 
latest reference being another example of data 
reanalysis). 

In spite of this weakness, we will assume that 
the model estimates of /' give a correct rep- 
resentation of ice volume anomalies around a 
tectonic-time- scale average. Ice- volume levels 
typical of the last interglacial then correspond 
to /' = —15 • 10^^ m^ in the model. The model 
prediction is an immediate but slow decrease in 
CO2 concentration (Figure[8| but without glacial 
inception before about 50,000 years (this is the 
Berger and Loutre prediction [H]!)- The particle 
filter also tells us that given the information at 
disposal (the model, the data, and the parameter 
priors) , it is not possible to provide a reliable es- 
timate of the evolution of climate beyond 50,000 
years. 

What about Ruddiman's hypothesis? Ruddi- 
man considers that humans perturbed climate's 
evolution around 8000 years ago. Therefore, we 
want to only consider data until that time, and 
see whether the model prediction differs to the 
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Ice volume anomaly {1 El 8 m3) 




C02 (ppmv) 



-20 -10 10 20 30 40 50 60 

Time (kyr) 

Figure 9: State estimate with the SM91 model, given 
data on CO2 and ice volume between 410 kyr BP and 
8 kyr BP (white) or lyr BP (grey). The subsequent 
prediction, with glacial inception in 50 kyr, is little 
affected by the data between 8 and kyr BP. This is 
opposed to Ruddiman's hypothesis 



parameter 


fixed value 


kfj, 


0.04 K / ppm /yr 


kg 


0.5 1/yr 


kR 


0.08 K / Wm^Vyr 


Ki 


l.e-4 1/yr 


Ke 


2.5e-4 yr"^ 



Table 1: Values of SM91 fixed parameters used both 
in the original publications and in the present article 



previous one. The experiment was carried out 
and the results are presented on Figure [9j The 
grey boxes provide the prediction with data as- 
similated until 8,000 years ago, and the white 
ones is the prediction with data assimilated until 
today. The two predictions are clearly undistin- 
guishable. Contrarily to Ruddiman, our model 
was therefore not "surprised" by the fact that 
CO2 continued to increase during the last 6,000 
years. 



5 Conclusion 

Behind this paper is the message that climate 
modelling is not and should not be a mere tech- 
nological question. Of course, general circula- 
tion models skillfully predict many complicated 
aspects of atmosphere and ocean dynamics; in 
that sense they are important and useful. Yet, 
they are but one aspect of the theoretical con- 
struct that underlies state-of-the-art knowledge 
of the climate system. Important questions are 
how we validate and calibrate climate models to 
provide the most informed predictions on climate 
change. 

Palaeo climates offer a premium playground to 
test the paradigms of complex system theory. We 
have been insistent on the fact that palaeocli- 
mate theory must rely on two pillars of modern 
applied mathematics: dynamical system theory 
and statistical decision theory. Along with the 
fact that palaeoclimate data have to be inter- 
preted and retrieved by skillful field scientists, 
their analysis turns to be a truly multidisci- 
plinary experience. The exceptionnally difficult 
challenges so posed are definitevely at the fron- 
tier of knowledge. 
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